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Abstract 

The FitzHugh-Nagumo equation has been investigated with a wide ar- 
ray of different methods in the last three decades. Recently a version of 
the equations with an applied current was analyzed by Champneys, Kirk, 
Knobloch, Oldeman and Sneyd IF using numerical continuation methods. 
They obtained a complicated bifurcation diagram in parameter space fea- 
turing a C-shaped curve of homoclinic bifurcations and a U-shaped curve 
of Hopf bifurcations. We use techniques from multiple time-scale dynam- 
ics to understand the structures of this bifurcation diagram based on geo- 
metric singular perturbation analysis of the FitzHugh-Nagumo equation. 
Numerical and analytical techniques show that if the ratio of the time- 
scales in the FitzHugh-Nagumo equation tends to zero, then our singular 
limit analysis correctly represents the observed CU-structure. Geometric 
insight from the analysis can even be used to compute bifurcation curves 
which are inaccessible via continuation methods. The results of our anal- 
ysis are summarized in a singular bifurcation diagram. 

1 Introduction 

1.1 Fast-Slow Systems 

Fast-slow systems of ordinary differential equations (ODEs) have the general 
form: 

dx 

ex = e—^f{x,y,€) (1) 

dy , ^ 
y = 9[x, y, e) 

where x £ M™, y G M" and < e ^ 1 represents the ratio of time scales. The 
functions / and g are assumed to be sufficiently smooth. In the singular limit 
e — >■ the vector field ([1]) becomes a differential-algebraic equation. The alge- 
braic constraint / = defines the critical manifold Cq = {{x,y) G x K" : 
f{x,y,0) = 0}. Where D^fip) is nonsingular, the implicit function theorem 
implies that there exists a map h{x) = y parametrizing Cq as a graph. This 
yields the impHcitly defined vector field y = g{h{y),y,0) on Co called the slow 
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flow. 



We can change ([T]) to the fast time scale t — rjt and let e — > to obtain the 
second possible singular limit system 

/(a^,2/,0) (2) 


We call the vector field ^ parametrized by the slow variables y the fast subsys- 
tem or the layer equations. The central idea of singular perturbation analysis 
is to use information about the fast subsystem and the slow flow to under- 
stand the full system ([ij. One of the main tools is Fenichel's Theorem (see 
[TBI [T71 [T51 [in])- It states that for every e sufflciently small and Co normally 
hyperbolic there exists a family of invariant manifolds for the flow ([IJ . The 
manifolds are at a distance 0(e) from Co and the flows on them converge to the 
slow flow on Co as e — > 0. Points p e Co where Dxf{p) is singular are referred 
to as fold pointfl 

Beyond Fenichel's Theorem many other techniques have been developed. 
More detailed introductions and results can be found in [12l [34] [2^ from a 
geometric viewpoint. Asymptotic methods are developed in [53] [53] whereas 
ideas from nonstandard analysis are introduced in [8]. While the theory is well 
developed for two-dimensional fast-slow systems, higher-dimensional fast-slow 
systems are an active area of current research. In the following we shall focus on 
the FitzHugh-Nagumo equation viewed as a three-dimensional fast-slow system. 



dx 
IE 
dy 

dt 



1.2 The FitzHugh-Nagumo Equation 

The FitzHugh-Nagumo equation is a simplification of the Hodgin-Huxley model 
for an electric potential of a nerve axon [30]. The first version was developed 
by FitzHugh |20l and is a two-dimensional system of ODEs: 

eii = V ■;r + u + P (3) 

V = — [V + — a) 
s 

A detailed summary of the bifurcations of ([3]) can be found in [44]. Nagumo et 
al. [13] studied a related equation that adds a diffusion term for the conduction 
process of action potentials along nerves: 

Ur = 5u^x + fa{u) -W+p 

Wr — e(u — jw) ^ 



^The projection of Co onto the x coordinates may liave more degenerate singularities than 
fold singularities at some of these points. 
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where fa{u) — u{u — a)(l — u) and p,j,5 and a are parameters. A good in- 
troduction to the derivation and problems associated with Q can be found in 
[2H] . Suppose we assume a travehng wave sohition to (|1]) and set u{x,t) = 
u{x + st) = u{t) and w{x,t) = w{x + st) = where s represents the wave 

speed. By the chain rule we get Ur = su', u^x — u" and Wr = sw' . Set v = u' 
and substitute into dH) to obtain the system: 

u' — V 

v' = ^{SV - fa{u) + W ~ p) (5) 

w' — -{u — jw) 
s 

System (O is the FitzHugh-Nagumo equation studied in this paper. Observe 
that a homoclinic orbit of (O corresponds to a traveling pulse solution of (g]). 
These solutions are of special importance in neuroscience |28| and have been 
analyzed using several different methods. For example, it has been proved that 
([S]) admits homoclinic orbits [221 13| for small wave speeds ( "slow waves" ) and 
large wave speeds ( "fast waves" ) . Fast waves are stable [33^ and slow waves are 
unstable [2T]. It has been shown that double-pulse homoclinic orbits [T5] are 
possible. If (O has two equilibrium points and heteroclinic connections exist, 
bifurcation from a twisted double heteroclinic connection implies the existence 
of multi-pulse traveling front and back waves [6] . These results are based on the 
assumption of certain parameter ranges for which we refer to the original papers. 
Geometric singular perturbation theory has been used successfully to analyze 
([5]) . In [32j the fast pulse is constructed using the exchange lemma [35j [31] [3] . 
The exchange lemma has also been used to prove the existence of a codimen- 
sion two connection between fast and slow waves in (s, e, a)-parameter space 
[55] . An extension of Fenichel's theorem and Melnikov's method can be em- 
ployed to prove the existence of heteroclinic connections for parameter regimes 
of (O with two fixed points [IS]. The general theory of relaxation oscillations 
in fast-slow systems applies to ([5]) (see e.g. [42l|26|) as does - at least partially 
- the theory of canards (see e.g. [46l l9l ITT] l39] ) . 

The equations ([5]) have been analyzed numerically by Champneys, Kirk, 
Knobloch, Oldeman and Sneyd [S] using the numerical bifurcation software 
AUTO [13l[T4]. They considered the following parameter values: 

7 = 1, " = ^' ^ = ^ 

We shall fix those values to allow comparison of our results with theirs. Hence 
we also write /i/io(u) — f{u). Changing from the fast time t to the slow time 
T and relabeling variables xi — u, X2 = v and y = w we get: 

exi — X2 

€±2 = ^{sx2 - xi{xi - - xi) + y - p) = l{sx2 - f{xi) + y - p){6) 

5 10 5 

y = -ixi-y) 
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From now on we refer to © as "the" FitzHugh-Nagumo equation. Investigating 
bifurcations in the (p, s) parameter space one finds C-shaped curves of homo- 
clinic orbits and a U-shaped curve of Hopf bifurcations; see Figure [T] Only 
part of the bifurcation diagram is shown in Figure [TJ There is another curve of 
homoclinic bifurcations on the right side of the U-shaped Hopf curve. Since ([SJ 
has the symmetry 

11 11 11 / 33 \ 

y-^-y, p^-^l-—j-p (7) 

we shall examine only the left side of the U-curve. The homoclinic C-curve is 
difficult to compute numerically by continuation methods using AUTO [T31 HI] 
or MatCont (35]. The computations seem infeasible for small values of e < 10"'^. 
Furthermore multi-pulse homoclinic orbits can exist very close to single pulse 
ones and distinguishing between them must necessarily encounter problems with 
numerical precision [5]. The Hopf curve and the bifurcations of limit cycles 
shown in Figure[l]have been computed using MatCont. The curve of homoclinic 
bifurcations has been computed by a new method to be described in Section [3.2l 




Figure 1: Bifurcation diagram of (|6]). Hopf bifurcations are shown in green, 
saddle-node of limit cycles (SNLC) are shown in blue and GH indicates a gen- 
eralized Hopf (or Bautin) bifurcation. The arrows indicate the side on which 
periodic orbits are generated at the Hopf bifurcation. The red curve shows 
(possible) homoclinic orbits; in fact, homoclinic orbits only exist to the left of 
the two black dots (see Section [3^ . Only part of the parameter space is shown 
because of the symmetry ([7]) . The homoclinic curve has been thickened to indi- 
cate that multipulse homoclinic orbits exist very close to single pulse ones (see 

Since the bifurcation structure shown in Figure[l]was also observed for other 
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excitable systems, Champneys et al. [S] introduced the term CU-system. Bifur- 
cation analysis from the viewpoint of geometric singular perturbation theory has 
been carried out for examples with one fast and two slow variables pTll^ l^HT] . 
Since the FitzHugh-Nagumo equation has one slow and two fast variables, the 
situation is quite different and new techniques have to be developed. Our main 
goal is to show that many features of the complicated 2-parameter bifurcation 
diagram shown in Figure[l]can be derived with a combination of techniques from 
singular perturbation theory, bifurcation theory and robust numerical methods. 
We accurately locate where the system has canards and determine the orbit 
structure of the homoclinic and periodic orbits associated to the C-shaped and 
U-shaped bifurcation curves, without computing the canards themselves. We 
demonstrate that the basic CU-structure of the system can be computed with 
elementary methods that do not use continuation methods based on collocation. 
The analysis of the slow and fast subsystems yields a "singular bifurcation dia- 
gram" to which the basic CU structure in Figure [T] converges as e — > 0. 

Remark: We have also investigated the termination mechanism of the C- 
shaped homoclinic curve described in [5]. Champneys et al. observed that the 
homoclinic curve does not reach the U-shaped Hopf curve but turns around and 
folds back close to itself. We compute accurate approximations of the homoclinic 
orbits for smaller values e than seems possible with AUTO in this region. One 
aspect of our analysis is a new algorithm for computing invariant slow manifolds 
of saddle type in the full system. This work will be described elsewhere. 



2 The Singular Limit 

The first step in our analysis is to investigate the slow and fast subsystems 
separately. Let e — > in ([6]); this yields two algebraic constraints that define 
the critical manifold: 

Co = |(a;i,a;2,y) e : X2 = y = Xi{xi - 1)(^ - xi)+p^ c{xi)^ 

Therefore Cq is a cubic curve in the coordinate plane X2 = 0. The parameter p 
moves the cubic up and down inside this plane. The critical points of the cubic 
are solutions of c'(xi) = and are given by: 

xi,± = (ll ± Vol) or numerically: « 0.6846, « 0.0487 

The points Xi^± are fold points with |c"(a;i ^)| 7^ since Cq is a cubic polynomial 
with distinct critical points. The fold points divide Cq into three segments 

Ci = {xi < a;i _}nCo, Cm = {xi- < xi < xi^+}nCo, Cr = {xi^+ < a;i}nCo 

We denote the associated slow manifolds by Ci_c, Cm,e and Cr.e- There are 
two possibilities to obtain the slow fiow. One way is to solve c{xi) = y for 



5 



Xi and substitute the result into the equation y = ^{xi — y). Alternatively 
differentiating y — c{xi) implicitly with respect to r yields y = xic'{xi) and 
therefore 

- y) = iic'(xi) ii = — - — -{xi-c{xi)) (8) 

s SC'(Xi) 

One can view this as a projection of the slow flow, which is constrained to the 
critical manifold in R^, onto the xi-axis. Observe that the slow flow is singu- 
lar at the fold points. Direct computation shows that the fixed point problem 
Xi = c{xi) has only a single real solution. This implies that the critical mani- 
fold intersects the diagonal y — xi only in a single point xl which is the unique 
equilibrium of the slow flow ([5]). Observe that q = (xj, 0,xj) is also the unique 
equilibrium of the full system ^ and depends on p. Increasing p moves the 
equilibrium from left to right on the critical manifold. The easiest practical 
way to determine the direction of the slow flow on Cq is to look at the sign of 
{xi — y). The situation is illustrated in Figure [2] 

Slow Flow in the plane x^-O 




_0_4^ \ \ \ \ \ \ \ 1 

-0.4 -0.2 0.2 0.4 0.6 0.8 1 1.2 



Figure 2: Sketch of the slow flow on the critical manifold Co 



2.1 The Slow Flow 

We are interested in the bifurcations of the slow flow depending on the parameter 
p. The bifurcations occur when x\ passes through the fold points. The values 
of p can simply be found by solving the equations c'{xi) = and c(xi) — Xi = 
simultaneously. The result is: 

p_ « 0.0511 and p+ « 0.5584 

where the subscripts indicate the fold point at which each equilibrium is located. 
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The singular time-rescaling r — sc'{xi)/t of the slow flow yields the desin- 
gularized slow flow 

— = xi - c{xi) = xi + — {xi- 1) (lOxi - 1) - p (9) 

Time is reversed by this rescaling on C; and Cr since s > and c'{xi) is 
negative on these branches. The desingularized slow flow ([9]) is smooth and has 
no bifurcations as p is varied. 



2.2 The Fast Subsystem 

The key component of the fast-slow analysis for the FitzHugh-Nagumo equation 
is the two-dimensional fast subsystem 

x'l = X2 

4 = \{sx2 - xi{xi - - xi) + y - p) (10) 

5 10 

where p > 0, s > are parameters and y is fixed. Since y and p have the same 
effect as bifurcation parameters we set p — y = p- We consider several fixed 
y-values and the effect of varying p (cf. Section 13. 2|) in each case. There are 
either one, two or three equilibrium points for ()10|) . Equilibrium points satisfy 
X2 = and lie on the critical manifold, i.e. we have to solve 

= Xl(xi-l)(^-.Tl)+p (11) 

for xi. We find that there are three equilibria for approximately = —0.1262 < 
p < 0.0024 = Pr, two equilibria on the boundary of this p interval and one 
equilibrium otherwise. The Jacobian of (|10p at an equilibrium is 

^(^^)=( ^(l-22"i + 30x?) I ) 

Direct calculation yields that for p ^ [phPr] the single equilibrium is a saddle. 
In the case of three equilibria, we have a source that lies between two saddles. 
Note that this also describes the stability of the three branches of the critical 
manifold C/, Cm and Cr- For s > the matrix A is singular of rank 1 if and only 
if 30x1 ~ 22a;i -f 1 = which occurs for the fold points xi^±. Hence the equilib- 
ria of the fast subsystem undergo a fold (or saddle-node) bifurcation once they 
approach the fold points of the critical manifold. This happens for parameter 
values pi and Pr- Note that by symmetry we can reduce to studying a single fold 
point. In the limit s = (corresponding to the case of a "standing wave") the 
saddle-node bifurcation point becomes more degenerate with A{xi) nilpotent. 
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Our next goal is to investigate global bifurcations of PH)) : we start with 
homoclinic orbits. For s = it is easy to see that (ITOl) is a Hamiltonian system: 



dH 
= t; — = X2 

0X2 



with Hamiltonian function 

H{x,,X2) - - ^ + ^ + — (1^^ 

We will use this Hamiltonian formulation later on to describe the geometry of 
homoclinic orbits for slow wave speeds. Assume that p is chosen so that (|12p 
has a homoclinic orbit XQ(t). We are interested in perturbations with s > 
and note that in this case the divergence of (ITOl) is s. Hence the vector field 
is area expanding everywhere. The homoclinic orbit breaks for s > and no 
periodic orbits are created. Note that this scenario does not apply to the full 
three-dimensional system as the equilibrium q has a pair of complex conjugate 
eigenvalues so that a Shilnikov scenario can occur. This illustrates that the 
singular limit can be used to help locate homoclinic orbits of the full system, 
but that some characteristics of these orbits change in the singular limit. 

We are interested next in finding curves in {p, s)-parameter space that rep- 
resent heteroclinic connections of the fast subsystem. The main motivation is 
the decomposition of trajectories in the full system into slow and fast segments. 
Concatenating fast heteroclinic segments and slow flow segments can yield ho- 
mochnic orbits of the full system [551I31I311[3S]. We describe a numerical strategy 
to detect heteroclinic connections in the fast subsystem and continue them in 
parameter space. Suppose that p G {pi,pr) so that pUj) has three hyperbolic 
equilibrium points xi, Xm and Xr- We denote by W^{xi) the unstable and by 
W^{xi) the stable manifold of xi. The same notation is also used for Xr and 
tangent spaces to Wi.) and W%) are denoted by T''(.) and T"(.). Recall that 
Xm is a source and shall not be of interest to us for now. Define the cross section 
Eby 

S = {(a;i,X2)eR2::ri = ^i±^}. 

We use forward integration of initial conditions in T" [xi ) and backward integra- 
tion of initial conditions in T'^{xr) to obtain trajectories 7''' and 7^ respectively. 
We calculate their intersection with E and define 

7/(p,s) 7+ n E, 7r(p, s) :=7"nE 

We compute the functions 7; and 7,- for different parameter values of (p, s) 
numerically. Heteroclinic connections occur at zeros of the function 

h{p,s) :=7/(p, s)-7r(p, s) 
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Once we find a parameter pair (po, so) such that h{pQ, sq) — 0, these parameters 
can be continued along a curve of heterochnic connections in (p, s) parameter 
space by solving the root-finding problem h{po + Si, so + 62) — for either 61 
or $2 fixed and small. We use this method later for different fixed values of y 
to compute heterochnic connections in the fast subsystem in (p, s) parameter 
space. The results of these computations are illustrated in Figure |31 There are 
two distinct branches in Figure [3] The branches are asymptotic to pi and pr 
and approximately form a "V" . From Figure [3] we conjecture that there exists 
a double heterochnic orbit for p « —0.0622. 



1.5 
1 

s 

0.5 





-0.15 -0.1 -0.05 0.05 



P 



Figure 3: Heterochnic connections for equation (jlOp in parameter space. 



Remarks: If we fix p = our initial change of variable becomes —y=p and 
our results for heterochnic connections are for the FitzHugh-Nagumo equation 
without an applied current. In this situation it has been shown that the hete- 
rochnic connections of the fast subsystem can be used to prove the existence of 
homoclinic orbits to the unique saddle equilibrium (0,0,0) (cf. |32]). Note that 
the existence of the heteroclinics in the fast subsystem was proved in a special 
case analytically [1 but Figure [3] is - to the best of our knowledge - the first ex- 
plicit computation of where fast subsystem heteroclinics are located. The paper 
|36j develops a method for finding heterochnic connections by the same basic 
approach we used, i.e. defining a codimension one hyperplane H that separates 
equilibrium points. 

Figure 13] suggests that there exists a double heterochnic connection for s — Q. 

Observe that the Hamiltonian in our case is H{xi,X2) = ^ + V{xi) where 
the function V{xi) is: 

PXI , ll{xif 



V{xi) 



5 100 150 20 



The solution curves of are given by X2 = ± a/2 (const. — V{xi)). The 
structure of the solution curves entails symmetry under refiection about the Xi- 
axis. Suppose p £ [pi,pr] and recall that we denoted the two saddle points of 
(jlOp by xi and Xr and that their location depends on p. Therefore, we conclude 
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that the two saddles xi and Xr must have a heterochnic connection if they he 
on the same energy level, i.e. they satisfy V{xi) — V{xr) = 0. This equation 
can be solved numerically to very high accuracy. 

Proposition 1. The fast subsystem of the FitzHugh-Nagumo equation for s = 
has a double heteroclinic connection for p = p* —0.0619259. Given a 
particular value y — yo there exists a double heteroclinic connection for p = 
P* + yo in the fast subsystem lying in the plane y = yo. 



2.3 Two Slow Variables, One Fast Variable 

From continuation of periodic orbits in the full system - to be described in 
Section 13.11 - we observe that near the U-shaped curve of Hopf bifurcations the 
a;2-coordinate is a faster variable than xi. In particular, the small periodic 
orbits generated in the Hopf bifurcation lie almost in the plane X2 = 0. Hence 
to analyze this region we set X2 = x^jt to transform the FitzHugh-Nagumo 
equation ([6]) into a system with 2 slow and 1 fast variable: 

X\ — X2 

e^X2 = i(sex2 - a;i(a;i - 1)(^ - Xl) + y -p) (14) 

y = -{xi~y) 
s 

Note that corresponds to the FitzHugh-Nagumo equation in the form (cf. 

my- 

Ur = He'^Uxx + f{u) -w +p , . 

Wr ~ e{u — w) 

Therefore the transformation X2 = 2^2 /e can be viewed as a rescaling of the 
diffusion strength by e^. We introduce a new independent small parameter 
6 = and then let 5 = — >■ 0. This assumes that 0(e) terms do not vanish in 
this limit, yielding the diffusion free system. Then the slow manifold Sq of (fTl]) 
is: 

5*0 = |(xi,X2,y) e R3 . ^ 1 -y+p)^ (16) 

Proposition 2. Following time rescaling by s, the slow flow of system (|14p on 
Sq in the variables {xi,y) is given by 

exi = f{xi)~y+p 
y = xi-y (17) 

In the variables {xi,X2) the vector field ()17p becomes 

Xl = X2 

ei2 = -'\ixi-f{xi)-p) + ^{f{xi)-e) (18) 
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Remark: The reduction to equations PT)) -(IT ^ suggests that is a three 
time-scale system. Note however that ([T4|) is not given in the three time-scale 
form (e^ii, ei2, is) = ihi{z),h2iz), h^iz)) for z = (zi, Z2, 2:3) e and /i^ : -)• 
R (i = 1,2,3). The time-scale separation in (IT71) - (IT51) results from the singular 
1/e dependence of the critical manifold Sq ; see pi 



Proof, (of Proposition [5]) Use the defining equation for the slow manifold ([T| 
and substitute it into xi ^ X2- A rescaling of time by t — > st under the 
assumption that s > yields the result p7)) . To derive (IT8| differentiate the 
defining equation of 5o with respect to time: 

^2 = — (i;i/'(a;i) - y) = — (x2/'(a;i) - y) 
se se 

The equations ij = -(xi — y) and y = — sea;2 -|- f{xi) +p yield the equations 
(UHl). ' □ 

Before we start with the analysis of (fT7| we note that detailed bifurcation 
calculations for P7|) exist. For example, Rocsoreanu et al. [H] give a detailed 
overview on the FitzHugh equation (|17p and collect many relevant references. 
Therefore we shall only state the relevant bifurcation results and focus on the 
fast-slow structure and canards. Equation (fT7|) has a critical manifold given by 
y — f{xi) + p — c{xi) which coincides with the critical manifold of the full 
FitzHugh-Nagumo system ([5]). Formally it is located in but we still denote 
it by Cq. Recall that the fold points are located at 

xi,± = (11 ± VoT) or numerically: « 0.6846, xi _ « 0.0487 

Also recall that the y-nullcline passes through the fold points at: 

p_ « 0.0511 and p+ « 0.5584 
We easily find that supercritical Hopf bifurcations are located at the values 



, , 2057 / 11728171 359e 509e2 

PH±(e) = ±\ \ (19) 

j-u.,±y , g^gQ y 182250000 1350 2700 27 ^ ' 

For the case e = 0.01 we get ^^,-(0.01) w 0.05632 and p/f,+(0.01) w 0.55316. 
The periodic orbits generated in the Hopf bifurcations exist for p e {ph, - , Ph. + ) ■ 
Observe also that ph,±{0) ~ P±\ so the Hopf bifurcations of P7|) coincide in 
the singular limit with the fold bifurcations in the one-dimensional slow fiow 
([H]). We are also interested in canards in the system and calculate a first order 
asymptotic expansion for the location of the maximal canard in (jl7p following 
[37] : recall that trajectories lying in the intersection of attracting and repelling 
slow manifolds are called maximal canards. We restrict to canards near the fold 
point {xi^-,c{xi.-)). 
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Proposition 3. Near the fold point (a^i,-, c(xi,_)) the maximal canard in (p, e) 
parameter space is given by: 

p(e) = xi,_~c(xi,_) + |e + 0(e3/2) 
o 

Proof. Let y = y ~ p and consider the shifts 

xi^xi+xi-, y y + c{xi-), p ^ p + xi^^ - c{xi-) 
to translate the equilibrium of (|17p to the origin when p — 0. This gives 



y IQ ^ij-y^ f{xi,y) 
y = <^{xi-y-p) = e{g{xi,y)~p) (20) 

Now apply Theorem 3.1 in 37 to find that the maximal canard of (|20p is given 
by: 

Shifting the parameter p back to the original coordinates yields the result. □ 

If we substitute e — 0.01 in the previous asymptotic result and neglect terms 
of order ©(e"^/^) then the maximal canard is predicted to occur for p fa 0.05731 
which is right after the first supercritical Hopf bifurcation at ph- ~ 0.05632. 
Therefore we expect that there exist canard orbits evolving along the middle 
branch of the critical manifold Cm.o.oi in the full FitzHugh-Nagumo equation. 
Maximal canards are part of a process generally referred to as canard explosion 
p!Ol [39l [7j- In this situation the small periodic orbits generated in the Hopf 
bifurcation at p = Ph,~ undergo a transition to relaxation oscillations within 
a very small interval in parameter space. A variational integral determines 
whether the canards are stable [39l [26] . 



Proposition 4. The canard cycles generated near the maximal canard point in 
parameter space for equation |J7p are stable. 

Proof. Consider the differential equation ([TTl) in its transformed form ([20|) . Ob- 
viously this will not affect the stability analysis of any limit cycles. Let xi{h) 
and Xm{h) denote the two smallest xi-coordinates of the intersection between 



Co := {(xi, y) e R2 : y = ^x\ - xf = ^(xi)} 

and the line y — h. Geometrically xi represents a point on the left branch and 
Xm a point on the middle branch of the critical manifold Cq. Theorem 3.4 in 
[5^ tells us that the canards are stable cycles if the function 

Jx,(h) clxi g[xi,(l){xi)) 
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is negative for all values h G (0, where xi = is the second fold point 

of Cq besides xi = 0. In our case we have 

Jxi{h) X-^xl+xl 

with xi{h) £ [—-^,0) and Xm{h) G (0, -^p]- Figure S] shows a numerical plot of 
the function R{h) for the relevant values of h which confirnis the required result. 
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Figure 4: Plot of the function R{h) for h G (0, (/)(^)]. 



Remark: We have computed an explicit algebraic expression for R'{h) with a 
computer algebra system. This expression yields R'{h) < for /i G (0, (/)(^^^)], 
confirming that R{h) is decreasing. 

□ 

As long as we stay on the critical manifold Cq of the full system, the analysis 
of the bifurcations and geometry of (|17p give good approximations to the dy- 
namics of the FitzHugh-Nagumo equation because the rescaling X2 = tx^ leaves 
the plane X2 = ^ invariant. Next we use the dynamics of the a;2-coordinate 
in system (|18p to obtain better insight into the dynamics when X2 7^ 0. The 
critical manifold Dq of ((T5| is: 

-Do = {(xi, X2) G : sx2d(x\) =x\- c[x{)\ 

We are interested in the geometry of the periodic orbits shown in Figure [5] that 
emerge from the Hopf bifurcation at Ph,-- Observe that the amplitude of the 
orbits in the xi direction is much larger that than in the a:2-direction. Therefore 
we predict only a single small excursion in the X2 direction for p slightly larger 
than Ph,- as shown in Figures 5(a) and 5(c) The wave speed changes the 



amplitude of this X2 excursion with a smaller wave speed implying a larger 
excursion. Hence equation (|17|) is expected to be a very good approximation 
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for periodic orbits in the FitzHugh-Nagumo equation with fast wave speeds. 
Furthermore the periodic orbits show two excursions in the relaxation regime 



after the canard explosion; see Figure 5(b) 



.0.01, p^O.OSS. 5^1.37 



, p^O.06, .37 



-D.005 
-0.D1 
-0.015 





.0.01 , p^o.osa, s^o.? 




(a) Small orbit near Hopf (b) Orbit after canard ex- (c) Different wave speed 
point (p = 0.058, s = 1.37) plosion (p = 0.06, s = 1.37) (p = 0.058, s = 0.2) 



Figure 5: Geometry of periodic orbits in the (xi, X2)-variables of the 2- variable 
slow subsystem (fT51) . Note that here X2 = is shown. Orbits have been 
obtained by direct forward integration for e = 0.01. 



3 The Full System 
3.1 Hopf Bifurcation 

The characteristic polynomial of the linearization of the FitzHugh-Nagumo 
equation ([6]) at its unique equilibrium point is 



^ ' 5s \ s / V 50 25 5 5 

Denoting P{X) — cq + ciA + C2A^ + csA^, a necessary condition for P to have 
pure imaginary roots is that cq = ciC2. The solutions of this equation can be 
expressed parametrically as a curve (j>{xl), s{xl)): 

s(x*? = 506(.^1) 
^ 1 + lOe - 22x3; + 30(a;*)2 

p{xl) = {xlf~l.l{xlf + l.l (21) 

Proposition 5. In the singular limit e — )■ the U-shaped bifurcation curves 
of the FitzHugh-Nagumo equation have vertical asymptotes given by the points 
P- » 0.0510636 andp^ w 0.558418 and a horizontal asymptote given by {{p, s) : 
p € and s = 0}. Note that at p± the equilibrium point passes through 

the two fold points. 

Proof. The expression for s{xl)'^ in (PT|) is positive when 1 + lOe — 22x1 + 
30{xlf < 0. For values of xl between the roots of 1 - 22a;i + 30{xl)'^ = 0, 
s(x*)2 — in (|2ip as e — 0. The values of p_ and p+ in the proposition are 
approximations to the value oip{xl) in pT|) at the roots of 1— 22x*-|-30(a;*)2 = 0. 
As e — 7> 0, solutions of the equation s{xl)'^ — c> in ([2T|) yield values of xl that 
tend to one of the two roots of 1 — 22.x^ + 30{xl)'^ = 0. The result follows. □ 



14 




Figure 6: Hopf bifurcation at p « 0.083, s = 1 and e = 0.01. The critical 
manifold Cq is shown in red and periodic orbits are shown in blue. Only the 
first and the last critical manifold for the continuation run are shown; not all 
periodic orbits obtained during the continuation are displayed. 

The analysis of the slow subsystems pT)) and ^TE\\ gives a conjecture about 
the shape of the periodic orbits in the FitzHugh-Nagumo equation. Consider 
the parameter regime close to a Hopf bifurcation point. From (|17p we expect 
one part of the small periodic orbits generated in the Hopf bifurcation to lie 
close to the slow manifolds Ci^e and Cm,e- Using the results about equation 
(|18p we anticipate the second part to consist of an excursion in the X2 direction 
whose length is governed by the wave speed s. Figure [6] shows a numerical con- 
tinuation in MatCont of the periodic orbits generated in a Hopf bifurcation 
and confirms the singular limit analysis for small amplitude orbits. 

Furthermore we observe from comparison of the xi and X2 coordinates of the 
periodic orbits in Figure |6(b)| that orbits tend to lie close to the plane defined 
by X2 = 0. More precisely, the X2 diameter of the periodic orbits is observed to 
be 0(e) in this case. This indicates that the rescaling of Section [^751 can help to 
describe the system close to the U-shaped Hopf curve. Note that it is difficult 
to check whether this observation of an 0(e)-diameter in the a:2-coordinate per- 
sists for values of e < 0.01 since numerical continuation of canard-type periodic 
orbits is difficult to use for smaller e. 



In contrast to this, it is easily possible to compute the U-shaped Hopf curve 
using numerical continuation for very small values of e. We have used this 
possibility to track two generalized Hopf bifurcation points in three parameters 
{p,s,e). The U-shaped Hopf curve has been computed by numerical continua- 
tion for a mesh of parameter values for e between 10~^ and 10~^ using MatCont 
[22] . The two generalized Hopf points GHf 2 on the left half of the U-curve were 
detected as codimension two points during each continuation run. The results 
of this "three-parameter continuation" are shown in Figure [T] 

The two generalized Hopf points depend on e and we find that their singular 



15 



£=0.01 



e=10" 



0.1706 0.1707 0.1 



0.1709 0.171 0.1711 0.1 



712 P 



(a) GHl 



£=10" 



£=0.01 



3 75 I , , , , , , , 

0.051 0.052 0.053 0.054 0.055 0.056 0.057 0.058 0.059 p 



(b) GH^ 



Figure 7: Tracking of two generalized Hopf points (GH) in (p, s, e)-parameter 
space. Each point in the figure corresponds to a diff'erent value of e. The point 
GHl in 7(a)| corresponds to the point shown as a square in Figure [T] and the 
point GH2 in |7(b)| is further up on the left branch of the U-curve and is not 
displayed in Figure [TJ 



limits in (p, s)-paranieter space are approximately: 

GH^ w (p = 0.171, s = 0) and GH^ ~ {p = 0.051, s = 3.927) 

We have not found a way to recover these special points from the fast-slow de- 
composition of the system. This suggests that codimension two bifurcations are 
generally diffcult to recover from the singular limit of fast-slow systems. 

Furthermore the Hopf bifurcations for the full system on the left half of the 
U-curve are subcritical between GiJ| and G-ffJ and supercritical otherwise. For 
the transformed system (fT4| with two slow and one fast variable we observed 
that in the singular limit ([T71) for — >■ the Hopf bifurcation is supercritical. 
In the case of e = 0.01 the periodic orbits for ([6]) and p7)) exist in overlapping 
regions for the parameter p between the p- values of GH'^ '^^ and GiJ2 °^- This 
result indicates that (|14l) can be used to describe periodic orbits that will interact 
with the homoclinic C-curve. 

3.2 Homoclinic Orbits 

In the following discussion we refer to "the" C-shaped curve of homoclinic bi- 
furcations of system ([S|) as the parameters yielding a "single-pulse" homoclinic 
orbit. The literature as described in Section [TT2] shows that close to single-pulse 
homoclinic orbits we can expect multi-pulse homoclinic orbits that return close 
to the equilibrium point multiple times. Since the separation of slow manifolds 
C.^e is exponentially small, homoclinic orbits of different types will always occur 
in exponentially thin bundles in parameter space. Values of e < 0.005 are small 
enough that the parameter region containing all the homoclinic orbits will be 
indistinguishable numerically from "the" C-curve that we locate. 
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The history of proofs of the existence of homochnic orbits in the FitzHugh- 
Nagumo equation is quite extensive. The main step in their construction is 
the existence of a "singular" homoclinic orbit 79. We consider the case when 
the fast subsystem has three equilibrium points which we denote by xi G C;, 
Xm € Cm and Xr S Cr- Rccall that xi coincides with the unique equilibrium 
q = (a;*,0,x*) of the full system for p < p^. A singular homoclinic orbit is 
always constructed by first following the unstable manifold of xi in the fast sub- 
system given hy y — xl- 




Figure 8: Homoclinic orbits as level curves of H{xi,X2) for equation (|12|) with 
y = a^i- 

First assume that s = 0. In this case the Hamiltonian structure - see Section 
I2.2l and equation - can be used to show the existence of a singular homoclinic 
orbit. Figure |5] shows level curves H{xi, X2) = H{x\, 0) for various values of p. 
The double heteroclinic connection can be calculated directly using Proposition 
[1] and solving x\-\- p* — p for p. 

Proposition 6. There exists a singular double heteroclinic connection in the 
FitzHugh-Nagumo equation for s = and p w —0.246016 = p* . 

Techniques developed in |45j show that the singular homoclinic orbits exist- 
ing for s — and p G {p*,p-) must persist for perturbations of small positive 
wave speed and sufficiently small e. These orbits are associated to the lower 
branch of the C-curve. The expected geometry of the orbits is indicated by 
their shape in the singular limit shown in Figure |8l The double heteroclinic 
connection is the boundary case between the upper and lower half of the C- 
curve. It remains to analyze the singular limit for the upper half. In this case, 
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a singular homoclinic orbit is again formed by following the unstable manifold 
of xi when it coincides with the equilibrium q = (a;^,0,a;*) but now we check 
whether it forms a heteroclinic orbit with the stable manifold of Xr- Then we 
follow the slow flow on Cr and return on a heteroclinic connection to C; for a 
different y-coordinate with y > xl and y < c{xi^+) = f{xi) + p. From there 
we connect back via the slow flow. Using the numerical method described in 
Section [2^ we first set y = xl; note that the location of q depends on the value 
of the parameter p. The task is to check when the system 

x[ = X2 

4 = + (22) 

has heteroclinic orbits from Ci to Cr with y = x^. The result of this compu- 
tation is shown in Figure [9] as the red curve. We have truncated the result at 
p = —0.01. In fact, the curve in Figure|9]can be extended to p = p* . Obviously 
we should view this curve as an approximation to the upper part of the C-curve. 

1.6 
1.6 
1.4 
1.3 
"3 1.2 
1.1 
1 

0.9 

0.8 

-0.02 0.02 0.04 0.06 0.08 

P 

Figure 9: Heteroclinic connections for equation (j22p in parameter space. The red 
curve indicates left-to-right connections for y = xl and the blue curves indicate 
right-to-left connections ioT y = xl + v with v ~ 0.125,0.12,0.115 (from top to 
bottom). 

If the connection from Cr back to Ci occurs with vertical coordinate xl+v, it 
is a trajectory of system (|22|) with y — xl+v. Figure |9] shows values of (p, s) at 
which these heteroclinic orbits exist for v ~ 0.125,0.12,0.115. An intersection 
between a red and a blue curve indicates a singular homoclinic orbit. Further 
computations show that increasing the value of v slowly beyond 0.125 yields 
intersections everywhere along the red curve in Figure [9] Thus the values of v 
on the homoclinic orbits are expected to grow as s increases along the upper 
branch of the C-curve. Since there cannot be any singular homoclinic orbits for 
p S {p- , p+ ) we have to find the intersection of the red curve in Figure IHl with 
the vertical line p = p-. Using the numerical method to detect heteroclinic 
connections gives: 
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Proposition 7. The singular homoclinic curve for positive wave speed termi- 
nates at p = and s « 1.50815 = s* on the right and at p = p* and s = on 
the left. 

In {p, s)-parameter space define the points: 

A=ip*,0), i? = (p_,0), (23) 

In Figure [10] we have computed the homochnic C-curve for values of e between 
10^^ and 5 • 10^^. Together with the singular limit analysis above, this yields 
strong numerical evidence for the following conjecture: 

Conjecture 1. The C- shaped homoclinic bifurcation curves converge to the 
union of the segments AB and AC as e — >■ 0. 

Remark 1: Figure 4 of Krupa, Sandstede and Szmolyan [38] shows a "wedge" 
that resembles shown in Figure (TUl The system that they study sets p = and 
varies a with a « 1/2. For a = 1/2 and p = 0, the equilibrium point q is 
located at the origin and the fast subsystem with y = has a double hetero- 
clinic connection at q to the saddle equilibrium (1,0,0) € Cr- The techniques 
developed in |38j use this double heteroclinic connection as a starting point. 
Generalizations of the results in [35] might provide a strategy to prove Conjec- 
ture [Tj rigorously, a possibility that we have not yet considered. However, we 
think that 1-homoclinic orbits in the regime we study come in pairs and that 
the surface of 1-homoclinic orbits in (p, s, e) space differs qualitatively from that 
described by Krupa, Sandstede and Szmolyan. 

Remark 2: We have investigated the termination or turning mechanism of 
the C-curve at its upper end. The termination points shown in Figure [Tj have 
been obtained by a different geometric method. It relies on the observation that, 
in addition to the two fast heteroclinic connections, we have to connect near C; 
back to the equilibrium point q to form a homoclinic orbit; the two heteroclinic 
connections might persist as intersections of suitable invariant manifolds but we 
also have to investigate how the flow near Ci^^ interacts with the stable man- 
ifold Wlq). These results will be reported elsewhere, but we note here that 
Pturn{e) P-. 

The numerical calculations of the C-curves for e < 10~^ are new. Numerical 
continuation using the boundary value methods implemented in AUTO |14j or 
MatCont [22] becomes very difficult for these small values of e [5]. Even com- 
puting with values e = O(10~^) using boundary value methods is a numerically 
challenging problem. The method we have used does not compute the homo- 
clinic orbits themselves while it locates the homoclinic C-curve accurately in 
parameter space. To motivate our approach consider Figure [TT] which shows the 
unstable manifold (q) for different values of s and fixed p. We observe that 
homoclinic orbits can only exist near two different wave speeds si and S2 which 
define the parameters where W"{q) C W^{Ci^e) or W"{q) C W'^{Cr,e)- Figure 
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Figure 10: Singular limit (e = 0) of the C-curve is shown in blue and parts of 
several C-curves for e > have been computed (red). 

[TT] displays how changes as s varies for the fixed value p = 0.05. If s 

differs from the points Si and S2 that define the lower and upper branches of 
the C-curvc for the given value of p, then increases rapidly on W'^{q) away 
from q. The changes in sign of xi on W"{q) identify values of s with homoclinic 
orbits. The two splitting points that mark these sign changes are visible in 
Figure [m Since trajectories close to the slow manifolds separate exponentially 
away from them, we are able to assess the sign of xi unambiguously on tra- 
jectories close to the slow manifold and find small intervals {p, si ± 10"^^) and 
{p, S2 ± 10""'^^) that contain the values of s for which there are homoclinic orbits. 

The geometry of the orbits along the upper branch of the C-curve is obtained 
by approximating it with two fast singular heteroclinic connections and parts of 
the slow manifolds Cr,e and Ciy, this process has been described several times in 
the literature when different methods were used to prove the existence of "fast 
waves" (see e.g. [2ill[32]). 

4 Conclusions 

Our results are summarized in the singular bifurcation diagram shown in Fig- 
ure 1121 This figure shows information obtained by a combination of fast-slow 
decompositions, classical dynamical systems techniques and robust numerical 
algorithms that work for very small values of e. It recovers and extends to 
smaller values of e the CU-structure described in |5] for the FitzHugh-Nagumo 
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(a) e = 0.01, p = 0.05, s e [0.1, 0.9] (b) e = 0.01, p = 0.05, s e [0.9, 1.5] 



Figure 11: Strong "splitting", marked by an arrow, of the unstable manifold 
W"{q) (red) used in the calculation of the homoclinic C-curve for small values 
of e. The critical manifold Co is shown in blue. The spacing in s is 0.05 for both 
figures. 

equation. The U-shaped Hopf curve was computed with an explicit formula, 
and the homoclinic C-curve was determined by locating transitions between dif- 
ferent dynamical behaviors separated by the homoclinic orbits. All the results 
shown as solid lines in Figure [T^] have been obtained by considering a singular 
limit. The lines AB and AC as well as the slow flow bifurcation follow from 
the singular limit e —> yielding the fast and slow subsystems of the FitzHugh- 
Nagumo equation ([6|). The analysis of canards and periodic orbits have been 
obtained from equations (fT7|l and ([T8|) where the singular limit was used 
(see Section [^T5|) . We have also shown the C- and U-curves in Figure [T^] as dot- 
ted lines to orient the reader how the results from Proposition [5] and Conjecture 
[T]fit in. 

We also observed that several dynamical phenomena are difficult to recover 
from the singular limit fast-slow decomposition. In particular, the codimension 
two generalized Hopf bifurcation does not seem to be observable from the singu- 
lar limit analysis. Furthermore the homoclinic orbits can be constructed from 
the singular limits but it cannot be determined directly from the fast and slow 
subsystems that they are of Shilnikov-type. 

The type of analysis pursued here seems to be very useful for other multiple 
time-scale problems involving multi-parameter bifurcation problems. In future 
work, we shall give a geometric analysis of the folding/turning mechanism of the 
homoclinic C-curve, a feature of this system we have not been able to determine 
directly from our singular limit analysis. That work relies upon new methods 
for calculating Cj^e and Cr.e which are invariant slow manifolds of "saddle-type" 
with both stable and unstable manifolds. 

We end with brief historical remarks. The references cited in this paper 
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Figure 12: Sketch of the singular bifurcation diagram for the FitzHugh-Nagumo equation (jH). The points A,B and C are 
defined in ([23]) . The part of the diagram obtained from equations (fT7 |) . ()18|) corresponds to the case "e^ = and e 7^ and 
small". In this scenario the canards to the right oi p = p- are stable (see Proposition 01). The phase portrait in the upper 
right for equation (flT)) shows the geometry of a small periodic orbit generated in the Hopf bifurcation of (fT7| . The two phase 
portraits below it show the geometry of these periodic orbits further away from the Hopf bifurcation for ([T7| . ([T8 |) . Excursions 
of the periodic orbits/canards for p > p- decrease for larger values of s. Note also that we have indicated as dotted lines the 
C-curve and the U-curve for positive e to allow a qualitative comparison with Figure [TJ 



discuss mathematical challenges posed by the FitzHugh-Nagumo equation, how 
these challenges have been analyzed and their relationship to general questions 
about multiple time-scale systems. Along the line AB in Figure[T2]we encounter 
a perturbation problem regarding the persistence of homoclinic orbits that can 
be solved using Fenichel theory |45| . The point A marks the connection be- 
tween fast and slow waves in (p, s)-parameter space which has been investigated 
in (e, s)-parameter space in [38]. We view this codimension 2 connectivity as one 
of the key features of the FitzHugh-Nagumo system. The perturbation problem 
for homoclinic orbits close to the line AC was solved using several methods and 
was put into the context of multiple time-scale systems in [211 131] , where the 
Exchange Lemma overcame difficulties in tracking W^{q) when it starts jump- 
ing away from Cr,e- This theory provides rigorous foundations that support our 
numerical computations and their interpretation. 
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